---
title: "IG Lobbying Analysis"
output: html_document
date: "2023-04-26"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(knitr)
library(readr)
library(tidyverse)
library(stargazer)
library(ggplot2)
library(gridExtra)
library(clubSandwich)
library(ggtext)
library(grid)
library(ggeffects)
library(lfe)
library(dotwhisker)
library(fect)
library(scales)
library(panelView)

#Not in function
'%!in%' <- function(x,y)!('%in%'(x,y))
```


```{r}
d <- read.csv("IG_Lobbying_Spending_Data_Markarian.csv")

```

```{r}
s <- d %>%
  pivot_longer(c(total_ig_cont_c, total_ig_cont_r))

ggplot(s, aes(x=control, y = value, fill=as.factor(name))) + 
  geom_boxplot() +
  labs(x = "Party Control", y = "Yearly Lobbying State Expenditures", title = "", fill="Position")  + 
  scale_x_discrete(labels = c("Democrats", "Divided", "Republican")) +
  theme_bw() +
  scale_fill_manual(values = c("grey90", "grey40"), labels = c("Gun Control IG", "Gun Rights IG"))
  
```

```{r}
d <- d %>%
  group_by(state_abb) %>%
  arrange(year) %>%
  mutate(total_fatalities = nonwhite_fatalities.tvp.b + white_fatalities.tvp.b,
         m_white_shooting = ifelse((white_fatalities.tvp.b >= total_fatalities*0.5) & total_fatalities>0, 1, 0),
         m_nonwhite_shooting = ifelse((nonwhite_fatalities.tvp.b >= total_fatalities*0.5) & total_fatalities>0, 1, 0),
         lagged.m_white_shooting =  dplyr::lag(m_white_shooting),
         lagged.m_nonwhite_shooting =  dplyr::lag(m_nonwhite_shooting),
         treated_white = ifelse(m_white_shooting==1 | lagged.m_white_shooting==1, 1, 0),
         treated_nonwhite = ifelse(m_nonwhite_shooting==1 | lagged.m_nonwhite_shooting==1, 1, 0),
         treated_white = ifelse(is.na(treated_white), 0, treated_white),
         treated_nonwhite = ifelse(is.na(treated_nonwhite), 0, treated_nonwhite),
         ig_c_2year = total_ig_cont_c + dplyr::lead(total_ig_cont_c),
         ig_r_2year = total_ig_cont_r + dplyr::lead(total_ig_cont_r)
         )
  
```

```{r}
m1 <- felm(ig_c_2year ~ white_fatalities.tvp.b + nonwhite_fatalities.tvp.b  | state_abb + year | 0 | state_abb  , d)
m2 <- felm(ig_r_2year ~ white_fatalities.tvp.b + nonwhite_fatalities.tvp.b  | state_abb + year | 0 | state_abb  , d)
m3 <- felm(ig_c_2year~ white_fatalities.tvp.b + nonwhite_fatalities.tvp.b + control + lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b | state_abb + year | 0 | state_abb  , d)
m4 <- felm(ig_r_2year ~ white_fatalities.tvp.b + nonwhite_fatalities.tvp.b + control + lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  | state_abb + year | 0 | state_abb  , d)
m5 <- felm(ig_c_2year ~ white_fatalities.tvp.b + nonwhite_fatalities.tvp.b + control + lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b + as.factor(state_abb)*as.numeric(year) | state_abb + year | 0 | state_abb  , d)
m6 <- felm(ig_r_2year~ white_fatalities.tvp.b + nonwhite_fatalities.tvp.b + control + lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b +  as.factor(state_abb)*as.numeric(year) | state_abb + year | 0 | state_abb  , d)
```

```{r}
stargazer(m1, m3, m5, covariate.labels = c("Number of White Fatalities", "Number of REM Fatalities",  "Divided Government" , "Republican Government", "Total Laws", "Neighbor White Fatalities", "Neighbor Nonwhite Fatalities", "Population Black", "Population Latino", "Population Other Race", "Population Median Income"),   dep.var.labels = c("Lobbying Spending in T + (T+1)") ,  omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect Pro-gun Control Interest Group State Lobbying Expenditures", type = "html", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table_IGcontrollob_apr2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize", digits = 0)
```

```{r}
stargazer(m2, m4, m6,  covariate.labels = c("Number of White Fatalities", "Number of REM Fatalities",  "Divided Government" , "Republican Government", "Total Laws", "Neighbor White Fatalities", "Neighbor Nonwhite Fatalities", "Population Black", "Population Latino", "Population Other Race", "Population Median Income"),     dep.var.labels = c("Lobbying Spending in T + (T+1)") ,  omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect Pro-gun Control Interest Group State Lobbying Expenditures", type = "html", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table_IGrightslob_apr2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize", digits = 0)
```

```{r}
#Plot
m1df <- broom::tidy(m1) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b")) %>%
  mutate(model = "Model 1")

m2df <- broom::tidy(m2) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b")) %>%
  mutate(model = "Model 2")

m3df <- broom::tidy(m3) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b")) %>%
  mutate(model = "Model 3")

m4df <- broom::tidy(m4) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b")) %>%
  mutate(model = "Model 4")

m5df <- broom::tidy(m5) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b")) %>%
  mutate(model = "Model 5")

m6df <- broom::tidy(m6) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b")) %>%
  mutate(model = "Model 6")


models1 <- rbind(m1df, m3df, m5df) %>%
  relabel_predictors("white_fatalities.tvp.b" = "Number of\nWhite Fatalities",
                     "nonwhite_fatalities.tvp.b" = "Number of\nREM Fatalities") 

models2 <- rbind(m2df, m4df, m6df) %>%
  relabel_predictors("white_fatalities.tvp.b" = "Number of\nWhite Fatalities",
                     "nonwhite_fatalities.tvp.b" = "Number of\nREM Fatalities")

models1$model <- factor(models1$model, levels = c("Model 1", "Model 3", "Model 5"))

models2$model <- factor(models2$model, levels = c("Model 2", "Model 4", "Model 6"))

g1 <- dwplot(models1, ci = 0.95,
       vline = geom_vline(
           xintercept = 0,
           colour = "grey60",
           linetype = 2
       ), dot_args = list(aes(shape = model))) +
    theme_bw(base_size = 11) + xlab("") + ylab("") +
    theme(
        plot.title = element_text(face = "bold"),
        legend.position = c(0.007, 0.01),
        legend.justification = c(0, 0),
        legend.background = element_rect(colour = "grey80"),
        legend.title.align = .5
    ) +
    labs(subtitle = "Gun Control Interest Group Lobbying\nSpending Over 2 Years") +
    scale_colour_grey(
        start = .1,
        end = .7,
        breaks = c("Model 1", "Model 3", "Model 5"),
        name = "Model",
        labels = c("TWFE", "TWFE & Controls", "TWFE, Controls, & Linear Trends")
    ) +
    scale_shape_discrete(
          breaks = c("Model 1", "Model 3", "Model 5"),
        name = "Model",
        labels = c("TWFE", "TWFE & Controls", "TWFE, Controls, & Linear Trends")
    ) +
    guides(
        shape = guide_legend(""), 
        colour = guide_legend("")) +
        theme(legend.direction="horizontal") 


g2 <- dwplot(models2, ci = 0.95,
       vline = geom_vline(
         xintercept = 0,
         colour = "grey60",
         linetype = 2
       ), dot_args = list(aes(shape = model))) +
  theme_bw(base_size = 11) + xlab("") + ylab("") +
  labs(subtitle = "Gun Rights Interest Group Lobbying\nSpending Over 2 Years") +
    scale_colour_grey(
        start = .1,
        end = .7,
        breaks = c("Model 2", "Model 4", "Model 6"),
        name = "Model",
        labels = c("TWFE", "TWFE & Controls", "TWFE, Controls, & Linear Trends")
    ) +
    scale_shape_discrete(
          breaks = c("Model 2", "Model 4", "Model 6"),
        name = "Model",
        labels = c("TWFE", "TWFE & Controls", "TWFE, Controls, & Linear Trends")
    ) +
    guides(
        shape = guide_legend(""), 
        colour = guide_legend("")) +
        theme(legend.position = "none", axis.text.y = element_blank()) 



g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

mylegend<-g_legend(g1)

g1 <- g1 + theme(legend.position="none")

lay = rbind(c(1,1,1,1,2,2,2),
            c(1,1,1,1,2,2,2),
            c(1,1,1,1,2,2,2),
            c(1,1,1,1,2,2,2),
            c(1,1,1,1,2,2,2),
            c(1,1,1,1,2,2,2),
            c(1,1,1,1,2,2,2),
            c(1,1,1,1,2,2,2),
            c(1,1,1,1,2,2,2),
            c(NA,NA,4,4,4,NA))

p3 <- grid.arrange(g1,g2, mylegend, layout_matrix=lay)

p3
```




# SI


```{r}
out.cfe <- fect(ig_c_2year  ~   m_white_shooting, 
     data = d, index = c("state_abb", "year"), 
method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, vartype = "jackknife")
out.cfe 
p1 <- plot(out.cfe, main = "", ylab = "ATT Gun Control Lobbying", xlab = "Time Since Majority White Shooting",  proportion =  .4)  +
  theme(axis.title = element_text(size = 10),
        axis.text.y = element_text(size = 9),
        axis.text.x = element_text(size = 9)) + 
  scale_y_continuous(labels = label_comma())


out.cfe <- fect(ig_r_2year  ~   m_white_shooting, 
     data = d, index = c("state_abb", "year"), 
method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, vartype = "jackknife")
out.cfe 
p2 <- plot(out.cfe, main = "", ylab = "ATT Gun Rights Lobbying", xlab = "Time Since Majority White Shooting",  proportion =  .4) +
  theme(axis.title = element_text(size = 10),
        axis.text.y = element_text(size = 9),
        axis.text.x = element_text(size = 9)) + 
  scale_y_continuous(labels = label_comma())


out.cfe <- fect(ig_c_2year  ~   m_nonwhite_shooting, 
     data = d, index = c("state_abb", "year"), 
method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, vartype = "jackknife")
out.cfe 
p3 <- plot(out.cfe, main = "", ylab = "ATT Gun Control Lobbying", xlab = "Time Since Majority REM Shooting",  proportion =  .6) +
  theme(axis.title = element_text(size = 10),
        axis.text.y = element_text(size = 9),
        axis.text.x = element_text(size = 9)) + 
  scale_y_continuous(labels = label_comma())


out.cfe <- fect(ig_r_2year  ~   m_nonwhite_shooting, 
     data = d, index = c("state_abb", "year"), 
method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, vartype = "jackknife")
out.cfe 
p4 <- plot(out.cfe, main = "", ylab = "ATT Gun Rights Lobbying", xlab = "Time Since Majority REM Shooting",  proportion =  .6) +
  theme(axis.title = element_text(size = 10),
        axis.text.y = element_text(size = 9),
        axis.text.x = element_text(size = 9)) + 
  scale_y_continuous(labels = label_comma())

grid.arrange(p1,p2,p3,p4)

#panelview(total_ig_cont_r  ~   treated_white, data = d, index = c("state_abb", "year"))

#panelview(total_ig_cont_r  ~   treated_nonwhite, data = d, index = c("state_abb", "year"))
```
# Unit specific linear time trends

```{r}
out.cfe <- fect(ig_c_2year  ~   m_white_shooting, 
     data = d, index = c("state_abb", "year"), 
method = "polynomial", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, degree = 1, vartype = "jackknife")
out.cfe 
p1 <- plot(out.cfe, main = "", ylab = "ATT Gun Control Lobbying", xlab = "Time Since Majority White Shooting",  proportion =  .4)  +
  theme(axis.title = element_text(size = 10),
        axis.text.y = element_text(size = 9),
        axis.text.x = element_text(size = 9)) + 
  scale_y_continuous(labels = label_comma())


out.cfe <- fect(ig_r_2year  ~   m_white_shooting, 
     data = d, index = c("state_abb", "year"), 
method = "polynomial", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, degree = 1,  vartype = "jackknife")
out.cfe 
p2 <- plot(out.cfe, main = "", ylab = "ATT Gun Rights Lobbying", xlab = "Time Since Majority White Shooting",  proportion =  .4) +
  theme(axis.title = element_text(size = 10),
        axis.text.y = element_text(size = 9),
        axis.text.x = element_text(size = 9)) + 
  scale_y_continuous(labels = label_comma())


out.cfe <- fect(ig_c_2year  ~   m_nonwhite_shooting, 
     data = d, index = c("state_abb", "year"), 
method = "polynomial", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, degree = 1, vartype = "jackknife")
out.cfe 
p3 <- plot(out.cfe, main = "", ylab = "ATT Gun Control Lobbying", xlab = "Time Since Majority REM Shooting",  proportion =  .6) +
  theme(axis.title = element_text(size = 10),
        axis.text.y = element_text(size = 9),
        axis.text.x = element_text(size = 9)) + 
  scale_y_continuous(labels = label_comma())


out.cfe <- fect(ig_r_2year  ~   m_nonwhite_shooting, 
     data = d, index = c("state_abb", "year"), 
method = "polynomial", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, degree = 1, vartype = "jackknife")
out.cfe 
p4 <- plot(out.cfe, main = "", ylab = "ATT Gun Rights Lobbying", xlab = "Time Since Majority REM Shooting",  proportion =  .6) +
  theme(axis.title = element_text(size = 10),
        axis.text.y = element_text(size = 9),
        axis.text.x = element_text(size = 9)) + 
  scale_y_continuous(labels = label_comma())

grid.arrange(p1,p2,p3,p4)

panelview(total_ig_cont_r  ~   m_white_shooting, data = d, index = c("state_abb", "year"), main = "Majority White Mass Shooting Treatment Status")

panelview(total_ig_cont_r  ~   m_nonwhite_shooting, data = d, index = c("state_abb", "year"), main = "Majority REM Mass Shooting Treatment Status")
```